home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / DASSL.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1996-07-24  |  9.7 KB  |  476 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include <cfloat>
  32. #include <cmath>
  33.  
  34. #include "DASSL.h"
  35. #include "f77-fcn.h"
  36. #include "lo-error.h"
  37.  
  38. extern "C"
  39. {
  40.   int F77_FCN (ddassl, DDASSL) (int (*)(const double&, double*,
  41.                     double*, double*, int&,
  42.                     double*, int*),
  43.                 const int&, double&, double*, double*,
  44.                 double&, const int*, const double&,
  45.                 const double&, int&, double*,
  46.                 const int&, int*, const int&,
  47.                 const double*, const int*, 
  48.                 int (*)(const double&, double*,
  49.                     double*, double*, const
  50.                     double&, double*, int*));
  51. }
  52.  
  53. static DAEFunc::DAERHSFunc user_fun;
  54. static DAEFunc::DAEJacFunc user_jac;
  55. static int nn;
  56.  
  57. DASSL::DASSL (void) : DAE ()
  58. {
  59.   stop_time_set = 0;
  60.   stop_time = 0.0;
  61.  
  62.   liw = 0;
  63.   lrw = 0;
  64.  
  65.   sanity_checked = 0;
  66.  
  67.   info.resize (15);
  68.  
  69.   for (int i = 0; i < 15; i++)
  70.     info.elem (i) = 0;
  71. }
  72.  
  73. DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
  74.   : DAE (state, time, f)
  75. {
  76.   n = size ();
  77.  
  78.   stop_time_set = 0;
  79.   stop_time = 0.0;
  80.  
  81.   liw = 20 + n;
  82.   lrw = 40 + 9*n + n*n;
  83.  
  84.   sanity_checked = 0;
  85.  
  86.   info.resize (15);
  87.  
  88.   for (int i = 0; i < 15; i++)
  89.     info.elem (i) = 0;
  90. }
  91.  
  92. DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv,
  93.       double time, DAEFunc& f)
  94.   : DAE (state, deriv, time, f)
  95. {
  96.   n = size ();
  97.  
  98.   stop_time_set = 0;
  99.   stop_time = 0.0;
  100.  
  101.   DAEFunc::set_function (f.function ());
  102.   DAEFunc::set_jacobian_function (f.jacobian_function ());
  103.  
  104.   liw = 20 + n;
  105.   lrw = 40 + 9*n + n*n;
  106.  
  107.   sanity_checked = 0;
  108.  
  109.   info.resize (15);
  110.  
  111.   for (int i = 0; i < 15; i++)
  112.     info.elem (i) = 0;
  113. }
  114.  
  115. void
  116. DASSL::force_restart (void)
  117. {
  118.   restart = 1;
  119.   integration_error = 0;
  120. }
  121.  
  122. void
  123. DASSL::set_stop_time (double t)
  124. {
  125.   stop_time_set = 1;
  126.   stop_time = t;
  127. }
  128.  
  129. void
  130. DASSL::clear_stop_time (void)
  131. {
  132.   stop_time_set = 0;
  133. }
  134.  
  135. int
  136. ddassl_f (const double& time, double *state, double *deriv,
  137.       double *delta, int& ires, double *, int *)
  138. {
  139.   ColumnVector tmp_deriv (nn);
  140.   ColumnVector tmp_state (nn);
  141.   ColumnVector tmp_delta (nn);
  142.  
  143.   for (int i = 0; i < nn; i++)
  144.     {
  145.       tmp_deriv.elem (i) = deriv [i];
  146.       tmp_state.elem (i) = state [i];
  147.     }
  148.  
  149.   tmp_delta = user_fun (tmp_state, tmp_deriv, time);
  150.  
  151.   if (tmp_delta.length () == 0)
  152.     ires = -2;
  153.   else
  154.     {
  155.       for (int i = 0; i < nn; i++)
  156.     delta [i] = tmp_delta.elem (i);
  157.     }
  158.  
  159.   return 0;
  160. }
  161.  
  162. int
  163. ddassl_j (const double& time, double *, double *, double *pd, const
  164.       double& cj, double *, int *)
  165. {
  166.   ColumnVector tmp_state (nn);
  167.   ColumnVector tmp_deriv (nn);
  168.  
  169.   // XXX FIXME XXX
  170.  
  171.   Matrix tmp_dfdxdot (nn, nn);
  172.   Matrix tmp_dfdx (nn, nn);
  173.  
  174.   DAEFunc::DAEJac tmp_jac;
  175.   tmp_jac.dfdxdot = &tmp_dfdxdot;
  176.   tmp_jac.dfdx    = &tmp_dfdx;
  177.  
  178.   tmp_jac = user_jac (tmp_state, tmp_deriv, time);
  179.  
  180.   // Fix up the matrix of partial derivatives for dassl.
  181.  
  182.   tmp_dfdx = tmp_dfdx +  cj * tmp_dfdxdot;
  183.  
  184.   for (int j = 0; j < nn; j++)
  185.     for (int i = 0; i < nn; i++)
  186.       pd [nn * j + i] = tmp_dfdx.elem (i, j);
  187.  
  188.   return 0;
  189. }
  190.  
  191. ColumnVector
  192. DASSL::do_integrate (double tout)
  193. {
  194.   ColumnVector retval;
  195.  
  196.   if (restart)
  197.     {
  198.       restart = 0;
  199.       info.elem (0) = 0;
  200.     }
  201.  
  202.   if (iwork.length () != liw)
  203.     iwork.resize (liw);
  204.  
  205.   if (rwork.length () != lrw)
  206.     rwork.resize (lrw);
  207.  
  208.   integration_error = 0;
  209.  
  210.   if (DAEFunc::jacobian_function ())
  211.     iwork.elem (4) = 1;
  212.   else
  213.     iwork.elem (4) = 0;
  214.  
  215.   double *px    = x.fortran_vec ();
  216.   double *pxdot = xdot.fortran_vec ();
  217.  
  218.   nn = n;
  219.   user_fun = DAEFunc::fun;
  220.   user_jac = DAEFunc::jac;
  221.  
  222.   if (! sanity_checked)
  223.     {
  224.       ColumnVector res = (*user_fun) (x, xdot, t);
  225.  
  226.       if (res.length () != x.length ())
  227.     {
  228.       (*current_liboctave_error_handler)
  229.         ("dassl: inconsistent sizes for state and residual vectors");
  230.  
  231.       integration_error = 1;
  232.       return retval;
  233.     }
  234.  
  235.       sanity_checked = 1;
  236.     }
  237.   
  238.   if (stop_time_set)
  239.     {
  240.       info.elem (3) = 1;
  241.       rwork.elem (0) = stop_time;
  242.     }
  243.   else
  244.     info.elem (3) = 0;
  245.  
  246.   double abs_tol = absolute_tolerance ();
  247.   double rel_tol = relative_tolerance ();
  248.  
  249.   if (initial_step_size () >= 0.0)
  250.     {
  251.       rwork.elem (2) = initial_step_size ();
  252.       info.elem (7) = 1;
  253.     }
  254.   else
  255.     info.elem (7) = 0;
  256.  
  257.   if (maximum_step_size () >= 0.0)
  258.     {
  259.       rwork.elem (2) = maximum_step_size ();
  260.       info.elem (6) = 1;
  261.     }
  262.   else
  263.     info.elem (6) = 0;
  264.  
  265.   double *dummy = 0;
  266.   int *idummy = 0;
  267.  
  268.   int *pinfo = info.fortran_vec ();
  269.   int *piwork = iwork.fortran_vec ();
  270.   double *prwork = rwork.fortran_vec ();
  271.  
  272. // again:
  273.  
  274.   F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo,
  275.                  rel_tol, abs_tol, idid, prwork, lrw,
  276.                  piwork, liw, dummy, idummy, ddassl_j));
  277.  
  278.   if (f77_exception_encountered)
  279.     (*current_liboctave_error_handler) ("unrecoverable error in dassl");
  280.   else
  281.     {
  282.       switch (idid)
  283.     {
  284.     case 1: // A step was successfully taken in intermediate-output
  285.             // mode. The code has not yet reached TOUT.
  286.     case 2: // The integration to TSTOP was successfully completed
  287.             // (T=TSTOP) by stepping exactly to TSTOP.
  288.     case 3: // The integration to TOUT was successfully completed
  289.             // (T=TOUT) by stepping past TOUT.  Y(*) is obtained by
  290.             // interpolation.  YPRIME(*) is obtained by interpolation.
  291.  
  292.       retval = x;
  293.       t = tout;
  294.       break;
  295.  
  296.     case -1: // A large amount of work has been expended.  (~500 steps).
  297.     case -2: // The error tolerances are too stringent.
  298.     case -3: // The local error test cannot be satisfied because you
  299.              // specified a zero component in ATOL and the
  300.          // corresponding computed solution component is zero.
  301.          // Thus, a pure relative error test is impossible for
  302.          // this component.
  303.     case -6: // DDASSL had repeated error test failures on the last
  304.          // attempted step.
  305.     case -7: // The corrector could not converge.
  306.     case -8: // The matrix of partial derivatives is singular.
  307.     case -9: // The corrector could not converge.  There were repeated
  308.          // error test failures in this step.
  309.     case -10: // The corrector could not converge because IRES was
  310.           // equal to minus one.
  311.     case -11: // IRES equal to -2 was encountered and control is being
  312.           // returned to the calling program.
  313.     case -12: // DDASSL failed to compute the initial YPRIME.
  314.     case -33: // The code has encountered trouble from which it cannot
  315.           // recover. A message is printed explaining the trouble
  316.           // and control is returned to the calling program. For
  317.           // example, this occurs when invalid input is detected.
  318.     default:
  319.       integration_error = 1;
  320.       break;
  321.     }
  322.     }
  323.  
  324.   return retval;
  325. }
  326.  
  327. Matrix
  328. DASSL::do_integrate (const ColumnVector& tout)
  329. {
  330.   Matrix dummy;
  331.   return integrate (tout, dummy);
  332. }
  333.  
  334. Matrix
  335. DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
  336. {
  337.   Matrix retval;
  338.   int n_out = tout.capacity ();
  339.  
  340.   if (n_out > 0 && n > 0)
  341.     {
  342.       retval.resize (n_out, n);
  343.       xdot_out.resize (n_out, n);
  344.  
  345.       for (int i = 0; i < n; i++)
  346.     {
  347.       retval.elem (0, i) = x.elem (i);
  348.       xdot_out.elem (0, i) = xdot.elem (i);
  349.     }
  350.  
  351.       for (int j = 1; j < n_out; j++)
  352.     {
  353.       ColumnVector x_next = do_integrate (tout.elem (j));
  354.  
  355.       if (integration_error)
  356.         return retval;
  357.  
  358.       for (int i = 0; i < n; i++)
  359.         {
  360.           retval.elem (j, i) = x_next.elem (i);
  361.           xdot_out.elem (j, i) = xdot.elem (i);
  362.         }
  363.     }
  364.     }
  365.  
  366.   return retval;
  367. }
  368.  
  369. Matrix
  370. DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out,
  371.           const ColumnVector& tcrit) 
  372. {
  373.   Matrix retval;
  374.   int n_out = tout.capacity ();
  375.  
  376.   if (n_out > 0 && n > 0)
  377.     {
  378.       retval.resize (n_out, n);
  379.       xdot_out.resize (n_out, n);
  380.  
  381.       for (int i = 0; i < n; i++)
  382.     {
  383.       retval.elem (0, i) = x.elem (i);
  384.       xdot_out.elem (0, i) = xdot.elem (i);
  385.     }
  386.  
  387.       int n_crit = tcrit.capacity ();
  388.  
  389.       if (n_crit > 0)
  390.     {
  391.       int i_crit = 0;
  392.       int i_out = 1;
  393.       double next_crit = tcrit.elem (0);
  394.       double next_out;
  395.       while (i_out < n_out)
  396.         {
  397.           int do_restart = 0;
  398.  
  399.           next_out = tout.elem (i_out);
  400.           if (i_crit < n_crit)
  401.         next_crit = tcrit.elem (i_crit);
  402.  
  403.           int save_output;
  404.           double t_out;
  405.  
  406.           if (next_crit == next_out)
  407.         {
  408.           set_stop_time (next_crit);
  409.           t_out = next_out;
  410.           save_output = 1;
  411.           i_out++;
  412.           i_crit++;
  413.           do_restart = 1;
  414.         }
  415.           else if (next_crit < next_out)
  416.         {
  417.           if (i_crit < n_crit)
  418.             {
  419.               set_stop_time (next_crit);
  420.               t_out = next_crit;
  421.               save_output = 0;
  422.               i_crit++;
  423.               do_restart = 1;
  424.             }
  425.           else
  426.             {
  427.               clear_stop_time ();
  428.               t_out = next_out;
  429.               save_output = 1;
  430.               i_out++;
  431.             }
  432.         }
  433.           else
  434.         {
  435.           set_stop_time (next_crit);
  436.           t_out = next_out;
  437.           save_output = 1;
  438.           i_out++;
  439.         }
  440.  
  441.           ColumnVector x_next = do_integrate (t_out);
  442.  
  443.           if (integration_error)
  444.         return retval;
  445.  
  446.           if (save_output)
  447.         {
  448.           for (int i = 0; i < n; i++)
  449.             {
  450.               retval.elem (i_out-1, i) = x_next.elem (i);
  451.               xdot_out.elem (i_out-1, i) = xdot.elem (i);
  452.             }
  453.         }
  454.  
  455.           if (do_restart)
  456.         force_restart ();
  457.         }
  458.     }
  459.       else
  460.     {
  461.       retval = integrate (tout, xdot_out);
  462.  
  463.       if (integration_error)
  464.         return retval;
  465.     }
  466.     }
  467.  
  468.   return retval;
  469. }
  470.  
  471. /*
  472. ;;; Local Variables: ***
  473. ;;; mode: C++ ***
  474. ;;; End: ***
  475. */
  476.